#define CBNTanalyser_cxx
#define PI 3.1415926
#include "CBNTanalyser.h"
#include <TH2.h>
#include <TStyle.h>
#include <string>
using namespace std;


void CBNTanalyser::Begin(TTree * tree/*tree*/)
{
   // The Begin() function is called at the start of the query.
   // When running with PROOF Begin() is only called on the client.
   // The tree argument is deprecated (on PROOF 0 is passed).
   TString option = GetOption();
   temproot = new TFile("amptAnalysis.root","RECREATE");
   look.open("look_logbook.txt");
   ptptpt = new TH1F("ptptpt","pt",1000,0,10);
   Ncol_init_low = new TH1F("Ncol_init_low","init pt < 1GeV",10,-0.5,9.5);
   Ncol_init_high = new TH1F("Ncol_init_high","init pt > 3GeV",10,-0.5,9.5);
   Ncol_final_low = new TH1F("Ncol_final_low","final pt < 1GeV",10,-0.5,9.5);
   Ncol_final_high = new TH1F("Ncol_final_high","final pt > 3GeV",10,-0.5,9.5);
   
   dpt_gain = new TH1F("dpt_gain","dpt_gain",120,-2,10);
   initial_phi_0 = new TH1F("initial_phi_0","Ncoll = 0",50,-3.1416,3.1416);
   initial_phi_1 = new TH1F("initial_phi_1","Ncoll = 1",50,-3.1416,3.1416);
   initial_phi_2 = new TH1F("initial_phi_2","Ncoll = 2",50,-3.1416,3.1416);
   initial_phi_3 = new TH1F("initial_phi_3","Ncoll >= 3",50,-3.1416,3.1416);

   f_initial_phi_0 = new TH1F("f_initial_phi_0","Ncoll = 0",50,-3.1416,3.1416);
   f_initial_phi_1 = new TH1F("f_initial_phi_1","Ncoll = 1",50,-3.1416,3.1416);
   f_initial_phi_2 = new TH1F("f_initial_phi_2","Ncoll = 2",50,-3.1416,3.1416);
   f_initial_phi_3 = new TH1F("f_initial_phi_3","Ncoll >= 3",50,-3.1416,3.1416);

   Final_init_pt = new TH2F("Final_init_pt","final pt vs initial pt",100,0,10,100,0,10);
   Final_init_phi = new TH2F("Final_init_phi","final vs initial phi",200,-3.1416,3.1416,200,-3.1416,3.1416);
   Final_init_phi_1 = new TH2F("Final_init_phi_1","Ncoll == 1",200,-3.1416,3.1416,200,-3.1416,3.1416);
   Final_init_phi_3 = new TH2F("Final_init_phi_3","Ncoll == 3",200,-3.1416,3.1416,200,-3.1416,3.1416);
   nVsDpt_1 = new TProfile("nVsDpt_1","final pt < .2GeV",20,-4,1);
   nVsDpt_2 = new TProfile("nVsDpt_2","final pt > 1GeV",20,-4,1);
   nVsDpt_3 = new TProfile("nVsDpt_3","final pt > 2GeV",20,-4,1);
   nVsDpt_4 = new TProfile("nVsDpt_4","final pt > 3GeV",20,-4,1);
   dptvsN_1 = new TProfile("dptvsN_1","final pt < .2GeV",10,-0.5,9.5);
   dptvsN_2 = new TProfile("dptvsN_2","final pt > 1GeV",10,-0.5,9.5);
   dptvsN_3 = new TProfile("dptvsN_3","final pt > 2GeV",10,-0.5,9.5);
   dptvsN_4 = new TProfile("dptvsN_4","final pt > 3GeV",10,-0.5,9.5);
   rptvsN_1 = new TProfile("rptvsN_1","final pt < .2GeV",10,-0.5,9.5);
   rptvsN_2 = new TProfile("rptvsN_2","final pt > 1GeV",10,-0.5,9.5);
   rptvsN_3 = new TProfile("rptvsN_3","final pt > 2GeV",10,-0.5,9.5);
   rptvsN_4 = new TProfile("rptvsN_4","final pt > 3GeV",10,-0.5,9.5);
   iptvsN_1 = new TProfile("iptvsN_1","initial pt < 0.2GeV",10,-0.5,9.5);
   iptvsN_2 = new TProfile("iptvsN_2","initial pt > 1GeV",10,-0.5,9.5);
   iptvsN_3 = new TProfile("iptvsN_3","initial pt > 2GeV",10,-0.5,9.5);
   iptvsN_4 = new TProfile("iptvsN_4","initial pt > 3GeV",10,-0.5,9.5);
   fptvsN_1 = new TH2F("fptvsN_1","initial pt < 0.2GeV",10,-0.5,9.5,200,0,10);
   initial_pt_t = new TH2F("initial_pt_t","initial_pt_t",200,0,20,200,0,10);
   pt_ncoll = new TH2F("pt_ncoll","special sample pt vs N",30,-0.5,29.5,100,0,10);

   dpt_max = new TH1F("dpt_max","max delta pt",600,-6,6);
   dpt_ratio = new TH1F("dpt_ratio","delta pt over final pt",200,-10,10);
   dpt_max_1 = new TH1F("dpt_max_1","max delta pt with low init pt",600,-6,6);
   dpt_ratio_ll = new TH1F("dpt_ratio_ll","delta pt over final pt with very low init pt",200,0,5);
   dpt_ratio_l = new TH1F("dpt_ratio_l","delta pt over final pt with low init pt",200,0,5);
   dpt_ratio_h = new TH1F("dpt_ratio_h","delta pt over final pt with high init pt",200,0,5);



   stringstream hname;
   stringstream hdescribe;
   stringstream profname1;
   stringstream profname2;
   stringstream profname3;
   stringstream profname4;
   stringstream profdescribe;
   float a[35];//pt
   float c[10];//# of collisiton
   float d[200];//cos(2*phi)
///////////////define pt bin from 0 to 6///////
   float ptbin[27];
   for(int i=0;i<20;i++)
	ptbin[i] = 0.1*i;
   for(int i=20;i<25;i++)
	ptbin[i] = 0.4 * i-6.;
   ptbin[25] = 4.;
   ptbin[26] = 6.; 
////////////////finish defining///////////////

   for(int i=0;i<20;i++)
	a[i] = 0.1*i;
   for(int i=20;i<30;i++)
	a[i] = 0.2*i-2.;
   a[30] = 4.;
   a[31] = 5.;
   a[32] = 6.;
   a[33] = 8.;
   a[34] = 10.;
   float b[201];
   for(int i=0;i<201;i++)
   {
   	b[i] = 3.1416/100.*i - 3.1416; 
	d[i] = 1./100*i-1.;
   }
   for(int i=0;i<11;i++)
	c[i] = i-0.5;
   float e[35];
   for(int i=0;i<35;i++)
	e[i] = 1./17*i-0.1;




   init_cos_init_pt_n = new TH3F("init_cos_init_pt_n","init v2 with initial pt",10,c,34,a,200,d);
   final_cos_init_pt_n = new TH3F("final_cos__init_pt_n","final v2 with initial pt",10,c,34,a,200,d);
   init_cos_final_pt_n = new TH3F("init_cos_final_pt_n","init v2 with final pt",10,c,34,a,200,d);
   final_cos_final_pt_n = new TH3F("final_cos_fianl_pt_n","final v2 with final pt",10,c,34,a,200,d);

   init_cos_init_pt_N = new TH3F("init_cos_init_pt_N","init v2 with initial pt",10,c,34,a,200,d);
   final_cos_init_pt_N = new TH3F("final_cos__init_pt_N","final v2 with initial pt",10,c,34,a,200,d);
   init_cos_final_pt_N = new TH3F("init_cos_final_pt_N","init v2 with final pt",10,c,34,a,200,d);
   final_cos_final_pt_N = new TH3F("final_cos_fianl_pt_N","final v2 with final pt",10,c,34,a,200,d);
   v2_n_i[0] = new TProfile("iv2_n_0","iv2_n_0",10,-0.5,9.5);
   v2_n_i[1] = new TProfile("iv2_n_1","iv2_n_1",10,-0.5,9.5);
   v2_n_i[2] = new TProfile("iv2_n_2","iv2_n_2",10,-0.5,9.5);
   v2_n_f[0] = new TProfile("fv2_n_0","fv2_n_0",10,-0.5,9.5);
   v2_n_f[1] = new TProfile("fv2_n_1","fv2_n_1",10,-0.5,9.5);
   v2_n_f[2] = new TProfile("fv2_n_2","fv2_n_2",10,-0.5,9.5);
   v2_n_Noptcut = new TProfile("v2_n_Noptcut","v2 vs n* without pt cut",10,-0.5,9.5);

   for(int i=0;i<6;i++)//i is to control the pt bin;
   {
	profname1.str("");
	profname1<<"v2_n_init_"<<i;
	profname2.str("");
	profname2<<"v2_n_final_"<<i;
	profdescribe.str("");
	profdescribe<<"v2 of "<<i;
	v2_N_initpt[i] = new TProfile(profname1.str().c_str(),profdescribe.str().c_str(),10,-0.5,9.5);
	v2_N_finalpt[i] = new TProfile(profname2.str().c_str(),profdescribe.str().c_str(),10,-0.5,9.5);
   }	
   for(int i=0;i<8;i++)//i is to control the pt bin;
   {
	profname3.str("");
	profname3<<"v2_pt_init_"<<i;
	profname4.str("");
	profname4<<"v2_pt_final_"<<i;
	profdescribe.str("");
	profdescribe<<"v2 of "<<i;
	v2_pt_initpt[i] = new TProfile(profname3.str().c_str(),profdescribe.str().c_str(),26,ptbin);
	v2_pt_finalpt[i] = new TProfile(profname4.str().c_str(),profdescribe.str().c_str(),26,ptbin);
   }	
		
   for(int i=0;i<11;i++)
   {
	hname.str("");
	hdescribe.str("");
	hname<<"final_phi_ipt_"<<i;
	hdescribe<<"When Ncoll = "<<i;
        final_phi_ipt[i] = new TH2F(hname.str().c_str(),hdescribe.str().c_str(),34,a,200,b);
   }
   for(int i=0;i<11;i++)
   {
	hname.str("");
	hdescribe.str("");
	hname<<"init_phi_ipt_"<<i;
	hdescribe<<"When Ncoll = "<<i;
        init_phi_ipt[i] = new TH2F(hname.str().c_str(),hdescribe.str().c_str(),34,a,200,b);
   }	
   for(int i=0;i<11;i++)
   {
	hname.str("");
	hdescribe.str("");
	hname<<"final_phi_fpt_"<<i;
	hdescribe<<"When Ncoll = "<<i;
        final_phi_fpt[i] = new TH2F(hname.str().c_str(),hdescribe.str().c_str(),34,a,200,b);
   }
   for(int i=0;i<11;i++)
   {
	hname.str("");
	hdescribe.str("");
	hname<<"init_phi_fpt_"<<i;
	hdescribe<<"When Ncoll = "<<i;
        init_phi_fpt[i] = new TH2F(hname.str().c_str(),hdescribe.str().c_str(),34,a,200,b);
   }

   for(int i=0;i<11;i++)
   {
	hname.str("");
	hdescribe.str("");
	hname<<"step_phi_pt_"<<i;
	hdescribe<<"When Ncoll = "<<i;
        step_phi_pt[i] = new TH2F(hname.str().c_str(),hdescribe.str().c_str(),34,a,200,b);
   }
  
   test_phi = new TH1F("test_phi","pt < 0.5",50,-3.1416,3.1416);
   test_phi_1 = new TH1F("test_phi_1","1 < pt < 2",50,-3.1416,3.1416);
   test_phi_2 = new TH1F("test_phi_2","pt > 2",50,-3.1416,3.1416);
   test_phi_3 = new TH1F("test_phi_3","0.5<pt<1",50,-3.1416,3.1416);
   itest_phi = new TH1F("itest_phi","pt < 0.5",50,-3.1416,3.1416);
   itest_phi_1 = new TH1F("itest_phi_1","1 < pt < 2",50,-3.1416,3.1416);
   itest_phi_2 = new TH1F("itest_phi_2","pt > 2",50,-3.1416,3.1416);
   itest_phi_3 = new TH1F("itest_phi_3","0.5<pt<1",50,-3.1416,3.1416);


   pt_tot = 0;
   n_tot = 0;

}

void CBNTanalyser::SlaveBegin(TTree *tree /*tree*/)
{
   // The SlaveBegin() function is called after the Begin() function.
   // When running with PROOF SlaveBegin() is called on each slave server.
   // The tree argument is deprecated (on PROOF 0 is passed).

   TString option = GetOption();

}

Bool_t CBNTanalyser::Process(Long64_t entry)
{
	fChain->GetTree()->GetEntry(entry);
	int size = mc_pt->size();
	double sumpt_init = 0;
	double sumpt_final = 0;	
	int sumnumber = 0;
	int sumncoll = 0;
	double sumcos_init = 0;
	double sumcos_final = 0;
	for(int i=0;i<size;i++)
	{
		if(mc_parent_index->at(i) < 0)
		{
////////////////////////parent things////////////////
			int Ncoll = 0;
			double pt_i = mc_pt->at(i);
			double eta_i = mc_eta->at(i);
			double phi_i = mc_phi->at(i);
			double d0_i = sqrt(pow(mc_vx_x->at(i),2)+pow(mc_vx_y->at(i),2));
			double z0_i = mc_vx_z->at(i);
			double px_i = mc_px->at(i);
			double py_i = mc_py->at(i);
			double t_i = mc_t->at(i);
			int child = i;
///////////////////////////Child things//////////////////
////////////////////////////////////
			if(pt_i<0.5)
				itest_phi->Fill(phi_i);
			if(pt_i<1 && pt_i>0.5)
				itest_phi_1->Fill(phi_i);
			if(pt_i<2 && pt_i>1)
				itest_phi_2->Fill(phi_i);
			if(pt_i<10 && pt_i>2)
				itest_phi_3->Fill(phi_i);
			while(mc_child_index->at(child) > 0)
			{
				child = mc_child_index->at(child);
				Ncoll++;
			}
			int id_i = mc_pdgId->at(i);
			int id_f = mc_pdgId->at(child);


			double pt_f = mc_pt->at(child);
			double px_f = mc_px->at(child);
			double py_f = mc_py->at(child);
			double eta_f = mc_eta->at(child);
			double phi_f = mc_phi->at(child);
			double d0_f = sqrt(pow(mc_vx_x->at(child),2)+pow(mc_vx_y->at(child),2));
			double z0_f = mc_vx_z->at(child);
//////////////////////////delta things//////////////
			ptptpt->Fill(pt_f);
//			if(pt_f<10)//3. && pt_i<1.0)
			{
				//look<<"the "<<i<<"th parton:"<<endl;
				//look<<"pt_i = "<<pt_i<<"\t pt_f = "<<pt_f<<endl;
				int ncol = 0;
				int dau = i;
				double pt_temp = mc_pt->at(i);
				double px_temp = 0;
                                double n_dpt=0;
                                double n_dpt_max = 0;
                                double phi_temp = mc_phi->at(i);
                                while(mc_child_index->at(dau) > 0)
                                {
                                        n_dpt = mc_pt->at(dau) - pt_temp;
                                        if(fabs(n_dpt) > fabs(n_dpt_max))
                                                n_dpt_max = n_dpt;
                                        pt_temp = mc_pt->at(dau);
                                        phi_temp = mc_phi->at(dau);
                                        dau = mc_child_index->at(dau);
//                                        step_phi_pt[Ncoll]->Fill(pt_temp,phi_temp);
					float cos2_temp = cos(2*phi_temp);
//					init_cos_init_pt_N->Fill(ncol,pt_temp,cos2_temp);
//					final_cos_init_pt_N->Fill(ncol,pt_temp,cos2_temp);
//					init_cos_final_pt_N->Fill(ncol,pt_temp,cos2_temp);
//					final_cos_final_pt_N->Fill(ncol,pt_temp,cos2_temp);
					v2_n_Noptcut->Fill(ncol,cos2_temp);
					if(pt_i<0.5)
					{
						v2_n_i[0]->Fill(ncol,cos2_temp);
						if(ncol == 0)
							test_phi->Fill(phi_temp);
					}
					if(pt_i<1 && pt_i>0.5)
						if(ncol == 0)
							test_phi_3->Fill(phi_temp);
						
					if(pt_i<2 && pt_i>1)
					{
						v2_n_i[1]->Fill(ncol,cos2_temp);
						if(ncol == 0)
							test_phi_1->Fill(phi_temp);
					}
					if(pt_i<10 && pt_i>2)
					{
						v2_n_i[2]->Fill(ncol,cos2_temp);
						if(ncol == 0)
							test_phi_2->Fill(phi_temp);
					}
					
					if(pt_f<0.5)
						v2_n_f[0]->Fill(ncol,cos2_temp);
					if(pt_f<2 && pt_f>1)
						v2_n_f[1]->Fill(ncol,cos2_temp);
					if(pt_f<10 && pt_f>2)
						v2_n_f[2]->Fill(ncol,cos2_temp);
					
                                        ncol++;
                                }
				{// Last collision cannot be recorded in while ()
				        n_dpt = mc_pt->at(dau) - pt_temp;
                                        if(fabs(n_dpt) > fabs(n_dpt_max))
                                                n_dpt_max = n_dpt;
                                        pt_temp = mc_pt->at(dau);
                                        phi_temp = mc_phi->at(dau);
					float cos2_temp = cos(2*phi_temp);
                                        step_phi_pt[Ncoll]->Fill(pt_temp,phi_temp);
					init_cos_init_pt_N->Fill(ncol,pt_temp,cos2_temp);
					final_cos_init_pt_N->Fill(ncol,pt_temp,cos2_temp);
					init_cos_final_pt_N->Fill(ncol,pt_temp,cos2_temp);
					final_cos_final_pt_N->Fill(ncol,pt_temp,cos2_temp);
					v2_n_Noptcut->Fill(ncol,cos2_temp);
					if(pt_i<0.5)
					{
						v2_n_i[0]->Fill(ncol,cos2_temp);
						if(ncol == 0)
							test_phi->Fill(phi_temp);
					}
					if(pt_i<2 && pt_i>1)
					{
						v2_n_i[1]->Fill(ncol,cos2_temp);
						if(ncol == 0)
							test_phi_1->Fill(phi_temp);
					}
					if(pt_i<10 && pt_i>2)
					{
						v2_n_i[2]->Fill(ncol,cos2_temp);
						if(ncol == 0)
							test_phi_2->Fill(phi_temp);
					}
					if(pt_i<1 && pt_i>0.5)
						if(ncol == 0)
							test_phi_3->Fill(phi_temp);

					if(pt_f<0.5)
						v2_n_f[0]->Fill(ncol,cos2_temp);
					if(pt_f<2 && pt_f>1)
						v2_n_f[1]->Fill(ncol,cos2_temp);
					if(pt_f<10 && pt_f>2)
						v2_n_f[2]->Fill(ncol,cos2_temp);
				}

				dpt_max_1->Fill(n_dpt_max);
				pt_ncoll->Fill(ncol,pt_temp);
			}
			double dpt = pt_f - pt_i;
/*
			if(Ncoll == 0 && dpt != 0)
				look<<"dpt = "<<dpt<<endl;
*/
			double rpt = 0;
			if(pt_i!=0)
				rpt = (pt_f-pt_i)/pt_i;
			double dphi = phi_f - phi_i;
			double ddphi = fabs(phi_f - phi_i);
			double iphicheck = atan2(py_i,px_i);
			if(ddphi > PI)
				dphi = dphi/ddphi*(2*PI - ddphi);
			double dphi2 = dphi * dphi;
//			if(pt_f>1.0 && Ncoll>0)
			if(1)//pt_i>1.5 && pt_f>1.5)
			{
				double rpt = 0;
				initial_pt_t->Fill(t_i,pt_i);
				if(pt_i<0.5)
					Ncol_init_low->Fill(Ncoll);
				if(pt_i>3.)
					Ncol_init_high->Fill(Ncoll);
				if(pt_f<0.5)
					Ncol_final_low->Fill(Ncoll);
				if(pt_f>3.)
					Ncol_final_high->Fill(Ncoll);
				Final_init_pt->Fill(pt_i,pt_f);
				Final_init_phi->Fill(phi_i,phi_f);
				if(Ncoll==1)
					Final_init_phi_1->Fill(phi_i,phi_f);
				if(Ncoll==3)
					Final_init_phi_3->Fill(phi_i,phi_f);
				if(pt_i<0.2)
				{
					rpt = dpt/(pt_i - 0.3087);
					rptvsN_1->Fill(Ncoll,rpt);
					dptvsN_1->Fill(Ncoll,dpt);
					iptvsN_1->Fill(Ncoll,pt_i);
					fptvsN_1->Fill(Ncoll,pt_f);
					nVsDpt_1->Fill(dpt,Ncoll);
				}
				if(pt_i>1)
				{
					rpt = dpt/(pt_i - 0.3087);
					rptvsN_2->Fill(Ncoll,rpt);
					dptvsN_2->Fill(Ncoll,dpt);
					iptvsN_2->Fill(Ncoll,pt_i);
					nVsDpt_2->Fill(dpt,Ncoll);
				}
				if(pt_i>2.)
				{
					rpt = dpt/(pt_i - 0.3087);
					rptvsN_3->Fill(Ncoll,rpt);
					dptvsN_3->Fill(Ncoll,dpt);
					iptvsN_3->Fill(Ncoll,pt_i);
					nVsDpt_3->Fill(dpt,Ncoll);
				}
				if(pt_i>3.)
				{
					rpt = dpt/(pt_i - 0.3087);
					rptvsN_4->Fill(Ncoll,rpt);
					dptvsN_4->Fill(Ncoll,dpt);
					iptvsN_4->Fill(Ncoll,pt_i);
					nVsDpt_4->Fill(dpt,Ncoll);
				}
				float cos2_init = cos(2*phi_i);
				float cos2_final = cos(2*phi_f);
			
				init_cos_init_pt_n->Fill(Ncoll,pt_i,cos2_init);
				final_cos_init_pt_n->Fill(Ncoll,pt_i,cos2_final);
				init_cos_final_pt_n->Fill(Ncoll,pt_f,cos2_init);
				final_cos_final_pt_n->Fill(Ncoll,pt_f,cos2_final);
				
/////////////////////////////////////////////////////////////////////////////
				if(pt_i>1. && Ncoll==6)
				{
					pt_tot += pt_f;
					n_tot += 1;
				}
				if(pt_f <0.5)
				{
					v2_N_finalpt[0]->Fill(Ncoll,cos2_init,1);
					v2_N_finalpt[1]->Fill(Ncoll,cos2_final,1);
				}

				if(1.<pt_f && pt_f<2.)
				{
					v2_N_finalpt[2]->Fill(Ncoll,cos2_init,1);
					v2_N_finalpt[3]->Fill(Ncoll,cos2_final,1);
				}
				if(pt_f>2.)
				{
					v2_N_finalpt[4]->Fill(Ncoll,cos2_init,1);
					v2_N_finalpt[5]->Fill(Ncoll,cos2_final,1);
				}

				if(pt_i<0.5)
				{
					v2_N_initpt[0]->Fill(Ncoll,cos2_init,1);
					v2_N_initpt[1]->Fill(Ncoll,cos2_final,1);
				}

				if(1.<pt_i && pt_i<2.)
				{
					v2_N_initpt[2]->Fill(Ncoll,cos2_init,1);
					v2_N_initpt[3]->Fill(Ncoll,cos2_final,1);
				}
				if(pt_i>2.)
				{
					v2_N_initpt[4]->Fill(Ncoll,cos2_init,1);
					v2_N_initpt[5]->Fill(Ncoll,cos2_final,1);
				}

				if(Ncoll == 0)
				{	
					v2_pt_initpt[0]->Fill(pt_i,cos2_init,1);
					v2_pt_initpt[1]->Fill(pt_i,cos2_final,1);
					v2_pt_finalpt[0]->Fill(pt_f,cos2_init,1);
					v2_pt_finalpt[1]->Fill(pt_f,cos2_final,1);
				}
				if(Ncoll == 1)
				{	
					v2_pt_initpt[2]->Fill(pt_i,cos2_init,1);
					v2_pt_initpt[3]->Fill(pt_i,cos2_final,1);
					v2_pt_finalpt[2]->Fill(pt_f,cos2_init,1);
					v2_pt_finalpt[3]->Fill(pt_f,cos2_final,1);
				}
				if(Ncoll == 2)
				{	
					v2_pt_initpt[4]->Fill(pt_i,cos2_init,1);
					v2_pt_initpt[5]->Fill(pt_i,cos2_final,1);
					v2_pt_finalpt[4]->Fill(pt_f,cos2_init,1);
					v2_pt_finalpt[5]->Fill(pt_f,cos2_final,1);
				}
				if(Ncoll > 2)
				{	
					v2_pt_initpt[6]->Fill(pt_i,cos2_init,1);
					v2_pt_initpt[7]->Fill(pt_i,cos2_final,1);
					v2_pt_finalpt[6]->Fill(pt_f,cos2_init,1);
					v2_pt_finalpt[7]->Fill(pt_f,cos2_final,1);
				}
////////////////////////////////////////////////////////////////////////////////////
				if(Ncoll <= 10)
				{
					final_phi_ipt[Ncoll]->Fill(pt_i,phi_f);
					init_phi_ipt[Ncoll]->Fill(pt_i,phi_i);
					final_phi_fpt[Ncoll]->Fill(pt_f,phi_f);
					init_phi_fpt[Ncoll]->Fill(pt_f,phi_i);
					int ncol = 0;
					int dau = i;
					double n_dpt=0;
					double n_dpt_max = 0;
					double pt_temp = mc_pt->at(i);
					double phi_temp = mc_phi->at(i);
					while(mc_child_index->at(dau) > 0)
                        		{
						n_dpt = mc_pt->at(dau) - pt_temp;
						dpt_ratio->Fill(n_dpt/pt_f);
						if(fabs(n_dpt) > fabs(n_dpt_max))
							n_dpt_max = n_dpt;
						pt_temp = mc_pt->at(dau);
						phi_temp = mc_phi->at(dau);
                                		dau = mc_child_index->at(dau);
						step_phi_pt[Ncoll]->Fill(pt_temp,phi_temp);
                                		ncol++;
                        		}

					{// Last collision cannot be recorded in while ()
					        n_dpt = mc_pt->at(dau) - pt_temp;
	                                        dpt_ratio->Fill(n_dpt/pt_f);
	                                        if(fabs(n_dpt) > fabs(n_dpt_max))
	                                                n_dpt_max = n_dpt;
	                                }
					dpt_max->Fill(n_dpt_max);
					if(Ncoll>1 && pt_i<0.2)
						dpt_ratio_ll->Fill(n_dpt_max/dpt);	
					if(Ncoll>1 && pt_i<0.5)
						dpt_ratio_l->Fill(n_dpt_max/dpt);	
					if(Ncoll>1 && pt_i>3.)
						dpt_ratio_h->Fill(n_dpt_max/dpt);	
					
				}
				if(pt_i>1. && pt_i<2.)
				{
					if(Ncoll==0)
						initial_phi_0->Fill(phi_i);
					if(Ncoll==1)
						initial_phi_1->Fill(phi_i);
					if(Ncoll==2)
						initial_phi_2->Fill(phi_i);
					if(Ncoll>2)
						initial_phi_3->Fill(phi_i);
				}
				if(pt_f>1. && pt_f<2.)
				{
					if(Ncoll==0)
						f_initial_phi_0->Fill(phi_i);
					if(Ncoll==1)
						f_initial_phi_1->Fill(phi_i);
					if(Ncoll==2)
						f_initial_phi_2->Fill(phi_i);
					if(Ncoll>2)
						f_initial_phi_3->Fill(phi_i);
				}

					
			}
	
		}			
	}
	

   return kTRUE;
}

void CBNTanalyser::SlaveTerminate()
{
   // The SlaveTerminate() function is called after all entries or objects
   // have been processed. When running with PROOF SlaveTerminate() is called
   // on each slave server.

}

void CBNTanalyser::Terminate()
{
   // The Terminate() function is the last function to be called during
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.
	look<<pt_tot<<endl<<n_tot<<endl;
        look.close();
        temproot->Write();
        temproot->Close();

}
